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ABSTRACT 

An algorithm for calculating acoustic intensities from a time-harmonic pressure field in an 
axisymmetric fluid region is presented. Acoustic pressures are computed in a mesh of NASTRAN 
triangular finite elements of revolution (TRIAAX) using an analogy between the scalar wave 
equation and elasticity equations. Acoustic intensities are then calculated from pressures and 
pressure derivatives taken over the mesh of TRIAAX elements. Intensities are displayed as vectors 
indicating the directions and magnitudes of energy flow at all mesh points in the acoustic field. A 
prolate spheroidal shell is modeled with axisymmetric shell elements (CONE AX) and submerged 
in a fluid region of TRIAAX elements. The model is analyzed to illustrate the acoustic intensity 
method and the usefulness of energy flow paths in the understanding of the response of 
fluid-structure interaction problems. The structural-acoustic analogy used is summarized for 
completeness. This study uncovered a NASTRAN limitation involving numerical precision issues 
in the CONEAX stiffness calculation causing large errors in the system matrices for nearly 
cylindrical cones. 

INTRODUCTION 

An acoustic intensity formulation for general, axisymmetric, fluid domains modeled with TRIAAX 
elements is presented. Numerical acoustic field solutions to fluid-structure interaction problems curt cntly yield 
acoustic pressure fields, which may be used to locate high acoustic pressure concentrations. The motivation for 
calculating and displaying acoustic intensities is to help visualize the energy flow paths which cause high pressure 
regions. The energy flow fields can then help to identify dominant power paths which flow between structure 
and fluid, and therefore the important radiating parts of a structure. 

The general problem of computing the interaction of an elastic structure with an acoustic fluid can be 
solved by combining a finite clement model of the structure with a fluid loading computed using boundary 
element [1-11], finite element [12-23], combined finite element/analytical [24-26], T-matrix [27-29], and 
approximate fluid loading [30-32] techniques. In the fluid finite element approach, the exterior fluid domain is 
modeled with finite elements truncated at a finite distance from the structure and terminated with an 
approximate radiation boundary condition to absorb outgoing waves. The principal computational trade-off 
between this approach and the boundary element approach is that the finite element approach yields large, 
banded matrices, whereas the boundary element approach yields smaller, densely-populated matrices. This 
trade-off sometimes favors the finite clement approach for long, slender structures like ships which are 
“naturally banded.” In addition, only the fluid finite element approach has directly available an explicit fluid 
mesh which can be used for graphical display of the wave motion through the fluid. Since a significant part of our 
interest involves the display of wave propagation through both structure and fluid, we therefore formulate this 
problem using the fluid finite element approach. The principal drawbacks to fluid finite element modeling are 
the need for an approximate radiation boundary condition at the outer fluid boundary, the requirements on mesh 
size and extent (sometimes leading to frequency-dependent fluid meshes [20]), and the difficulty of generating 
the fluid mesh. 
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Direct frequency response N ASTRAN [33] solutions for axisymmetric regions are described in general 
and for domains that are defined using both structural and fluid subregions. Structural regionsare modeled us g 
sta^d^rd^ Ructuraf f irhte elements (CONEAX, TRIAAX, TRAPAX); fluid regions are modeled with triangular 
elements of revolution (TRIAAX) using an analogy relating the Helmholtz equation to the elasticity equa ions 

is described id dejad b, Eveline 134 and d • — ” 
completeness here. The modeling of fluid-structure interaction between fluid and structure domains 
defined, as well as the application of acoustic boundary conditions to fluid models. 

The acoustic intensity formulation includes the definition of the intensity quantity, the ^Igonthm used 

to calculate fluid particle velocities using the pressure distribution in a f n f ral , ^SefomuZn kafbeen 
calculation of the acoustic intensity vector quantity from pressures and velocities. The formulation has been 
implemented in the program AcINT (Acoustic INTensity), which functions as a post-processor to NASIttAN. 
Aifexample a submerged prolate-spheroidal shell with two sets of boundary conditions, is analyzed for a ring 
load The resulting acoustic intensity fields are displayed for a given excitation frequency to illustrate the energy 
flow paths which result. The acoustic intensity vector plots show the utility of the method in identifying dom 
power paths in fluid-structure interaction problems. 

STRUCTURAL-ACOUSTIC ANALOGY 

From an engineering point of view, it is convenient to be able to make use of existing general purpose 
finite element codes for analyzing structural-acoustic problems. Finite element codes are widely available 
versatile reliable, well supported, and an abundance of pre- and postprocessors may be used with t em. us 
ZL in .his sectJTanalogy [34] between , he equa, ionsof 

This analogy allows the coupled structural acoustic problem to be solved with standard finite elemen 
NASTRAN. 


The pressure p in an acoustic field satisfies the wave equation 

V>P = 3> 


( 1 ) 


where V 2 is the Laplacian operator, p is the dynamic fluid pressure, c is the wave speed, and dots denote partial 
differentiation with respect to time. 


On the other hand, the x-component of the Navier equations of elasticity, which are the equations 
solved by structural analysis computer programs, is 
A + 2G 
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G - G 

where u v, and w are the Cartesian components of displacement, k is a Lame elastic constant, G is the shear 

modulus, f x is the x-component of body force per unit volume (e.g., gravity), q is the mass density, and commas 

denote partial differentiation. 

A comparison of Eqs. 1 and 2 indicates that elastic finite elements can be used to model scalar pressure 

fields if we let u, the x-component of displacement, represent p, set v = w = 0 everywhere, x - , an j. 

For three-dimensional analysis, the engineering constants consistent with this last requirement are 

E e = 10 20 G e , pb = (3) 

where the element shear modulus G e can be selected arbitrarily. The subscript “e” has been added to these 
constants to emphasize that they are merely numbers assigned to the elements. 

A variety of boundaiy conditions may also be imposed. At a pressure-release boundary, p = 0 is 
enforced explicitly like other displacement boundary conditions. For gradient conditions, the pressure gradien 
dp/dn is enforced at a boundary point by applying a “force” to the unconstrained DOF at that point equal to 
G e A6p/6n , where A is the area assigned to the point and n is the outward normal from the fluid region. For 
example, the plane wave absorbing boundary condition 

ap = _p (4) 

an c 
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is enforced by applying to each point on the outer fluid boundary a “force” given by (G e A/c)p - Since this 
“force” is proportional to the first time derivative of the fundamental solution variable p, this boundary Condition 
isimposed in the analogy by attaching to the fluid DOF a “dashpot” of constant G c A/c . The Neumann condition 
dp/dn = 0 is the natural boundary condition under this analogy. The next higher order local radiation 
boundary condition, the curved wave absorbing boundary condition [23,36] 

dp £ £ ( 5 ) 

an c r ’ 

where r is the radius of the boundary, is enforced under the analogy by attaching in parallel both a “dashpot” and a 
“spring” between each boundary point and ground. 

At a fluid-structure interface (an accelerating boundary), momentum and continuity considerations 
require that 

ap (6) 

-f- — eu n , 
an 

where n is the normal at the interface, q is the mass density of the fluid, and u n is the normal component of fluid 
particle acceleration. Under the analogy, this condition is enforced by applying to the fluid DOF at a 
fluid-structure interface a “force” given by - (G t pA)u„ . 


Tb summarize, the wave equation, Eq. 1, can be solved with elastic finite elements if the 
three-dimensional region is modeled with 3-D solid finite elements having material properties given by Eq. 3, 
and only one of the three Cartesian components of displacement is retained to represent the scalar variable p. In 
Cartesian coordinates, any of the three components can be used. The solution of axisymmetnc problems in 
cylindrical coordinates follows the same approach except that the z-component ol displacement is the on y one 
which can be used to represent p. 


FINITE ELEMENT FORMULATION OF FLUID-STRUCTURE INTERACTION 

There are two fundamental fluid-structure interaction problems of interest in structural acoustics: 
acoustic radiation, in which a submerged elastic body is subjected to a mechanical excitation applied to the 
structure, and scattering, in which the structure is subjected to an external incident pressure loading. For general 
time-dependent problems, the excitation is an arbitrary function of time, whereas in the time-harmonic case ol 
interest here, the excitation has a single circular frequency a» . 

Although our specific interest here is the time-harmonic case, we summarize the theory [22,35] for the 
more general case, which includes an incident loading as well. The radiation problem will be covered as a special 
case. For scattering, we assume, without loss of generality, that the incident wave propagates in the negative z 
direction. The speed of such propagation is c, the speed of sound in the fluid. 

Within the fluid region, the total pressure p satisfies the wave equation, Eq. 1. Since the incident 
free-field pressure p, is known, it is convenient to decompose the total pressure p into the sum of incident and 

scattered pressures 

(7) 

p = Pt + Ps, K } 

each of which satisfies the wave equation. 

We now formulate the problem for finite element solution. Consider an arbitrary, submerged, 
three-dimensional elastic structure subjected to either internal time-dependent loads or an external 
time-dependent incident pressure. If the structure is modeled with finite elements, the resulting matrix equation 
of motion for the structural degrees of freedom (DOF) is 

Mu + Bu + Ku = F - GAp, 

where M, B, and K are the structural mass, viscous damping, and stiffness matrices (dimension s x s), respectively, 
u is the displacement vector for all structural DOF (wet and dry) in terms of the coordinate systems selected by 
the user (s x r), F is the vector of applied mechanical forces applied to the structure (s x r), G is the rectangular 
transformation matrix of direction cosines to transform a vector of outward normal forces at the wet points to a 
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vector of forces at all points in the coordinate systems selected by the user (s x t). A is the diagonal area matrix foi 
tS wet surface (f x f),p is the vector of total fluid pressures (incident + scattered) applied at the wet grid | o mis (1 
^nd Its denote' different, at, on w„h respect ,0 time. The pressure p is a-med I = in ^mp = In 

the above dimensions, s denotes the total number of independent structural DOF (wet rnfirs, order finite 
nnmher of fluid DOF (the number of wet points), and r denotes the number of load cases. It first order 
ekmmttsareused for the suXe discretisation, surface areas, normals, and the transformation mat™ G can be 
obtained from the calculation of the load vector resulting from an outwardly directed static unit pressure load on 
the structure’s wet surface. The matrix product GA can then be interpreted as the matrix which converts a vccto 
of negative fluid pressures to structural loads in the global coordinate system. The last two equations can be 

combined to yield 

Mu + Bu + Ku 4- GAp s “ F GAp,. 


A finite element model of the fluid region (with scattered pressure p. as the unknown) results in a 
matrix equation of the form 

Qp. + Cp 5 + Hp s = F*\ (10) 

where p, is the vector of scattered fluid pressures at the grid points of the fluid region, Q and H are the fluid 
“inertia” and “stiffness” matrices (analogous to M and K for structures), C is the “damping” matrix arising from 
the radiation boundary condition (Eq. 4), and P p) is the “loading” applied to fluid DOF due to the fluid-structure 
interface condition, Eq. 6. Using the analogy described in the preceding section, elastic finite elements ca 
used to model both structural and fluid regions. Material constants assigned to the e astic elements used to 
model the fluid .are specified according to Eq. 3. In three dimensions, elastic sohddements areusede.g 
isoparametric bricks (IHEXi) for general 3-D analysis or solids of revolution (TRIAAX, 1RAPAX) lor 
axisymmetric analysis). 


At the fluid-structure interface, Eqs. 6 and 7 can be combined to yield 


— = e(u ni - u n ), 
dn 


( 11 ) 


where n is the outward unit normal, and u„, and u„ arc, respectively, the incident and total outward normal 
components of fluid particle acceleration at the interlace. Thus, from the analogy, we impose the fluid-structure 
interface condition by applying a “load” to each interface fluid point given by 

F® = -<?G e A(u ni - u„), ' 

where the first minus sign is introduced since, in the coupled problem, we choose n as the outward normal from 
the structure into the fluid, making n an inward normal for the fluid region. The normal displacements u n a 
related to the total displacements u by the same rectangular transformation matrix G used above: 

u„ = G T u, (13) 

where the superscript T denotes the matrix transpose. Eqs. 10, 12, and 13 can be combined to yield 

Qp. + Cp s + Hp. - pG e (GA) T u = - pGcAiin,. ( ) 

Since the fluid-structure coupling terms in Eqs. 9 and 14 are nonsymmetric, we symmetrize the problem [21] by 
using a new fluid unknown q such that 

(15) 


l 

= | Psdt, 


q = Ps 


If Eq. 14 is integrated in time, and the fluid element “shear modulus” C, is chosen as 
-1 

G - ' T' 

the overall matrix system describing the coupled problem can be written as 


(16) 
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( 17 ) 


M O' 


B (GA)' 

H 


'k o' 

ful 

_ (F-(GA)p/ 

0 0 

4 

(GA) 1 C 

k 

■ ~\- 

0 H 

k 

1 Av "' 


where v ni (= u ni ) is the outward normal component of incident fluid particle velocity. 


The new variable C{ is, except lor a multiplicative constant, the velocity potential 0 , since 

p = ~p0. 

Eq. 17 could also be recast in terms of 0 rather than q as the fundamental fluid unknown, but no particular 
advantage would result. In fact, the use ol cj rather than 0 has the practical advantage that the fluid pressure can 
be recovered directly from the finite element program as the time derivative (velocity) of the unknown q. 

lo summarize, both structural and fluid regions are modeled with finite elements. For the fluid region, 
the material constants assigned to the finite elements are 

- 10 20 - 1 - 1 

E e « , G c = — , v e = unspecified, g c = — y, (19) 

q q Q c 

where E e , G c , v c , and p c are the Young’s modulus, shear modulus, Poisson’s ratio, and mass density, 
respectively, assigned to the fluid finite elements. The properties q and c above are the actual density and sound 

speed for the fluid medium. The radiation boundary condition used is the plane wave approximation, Eq. 4, which 
appears to be adequate if the outer fluid boundary is sufficiently far from the structure [20]. With this boundary 
condition, matrix C in Eq. 17 arises from dashpots applied at the outer fluid boundary with damping constant 
- A/(qc) at each grid point to which the area A has been assigned. At the fluid-structure interface, matrix G A is 

entered using the areas (or areal direction cosines) assigned to each wet degree of freedom. (Recall that GA can 
be interpreted as the matrix which converts a vector of negative fluid pressures to structural loads in the global 
coordinate system.) 

For radiation problems, the right-hand side of Eq. 17 can be simplified further since the incident 
pressure p, is zero, and we obtain 
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We note that the structural and fluid unknowns are not sequenced as perhaps implied by the partitioned 
form of Eq. 20. The coupling matrix G A is quite sparse and has nonzeros only for matrix rows associated with the 
structural DOF at the fluid-structure interface and columns associated with the coincident fluid points. Thus, the 
grid points should be sequenced for minimum matrix bandwidth or profile as if the structural and fluid meshes 
comprised a single large mesh. As a result, the structural and fluid grid points will, in general, be interspersed in 
their numbering, and the system matrices will be sparse and banded. 

ACOUSTIC INTENSITY CALCULATIONS 

The procedure for solving for the acoustic intensity field in an axisymmetric fluid finite element model 
using NASTRAN [33] and the acoustic intensity post-processor AcINT is: 

• run NASTRAN on a dynamically loaded finite element model of structural elements 
and fluid elements using direct frequency response analysis, and generate resultant 
nodal pressures for the fluid region(s); 

• run AcINT using the output from NASTRAN to calculate nodal fluid velocities and 
acoustic intensities. 

The nodal pressures, actually the cosine coefficients of nodal pressures, are computed by NASTRAN 
in response to the AXISYMMETRIC = COSINE command in the case control part of the input data. Since only 
cosine coefficients are requested, the 2, 4, and 6 DOF are removed by NASTRAN for harmonic zero. The 1 
and 5 DOF must be constrained by the user, leaving the 3 DOF to represent the scalar fluid velocity potential. 
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The velocity potential integrated over time is the fluid pressure. Therefore, in the analogy, nodal pressures are 
obtained bythe case control command VELOCITY = ALL (or SETID). The resultant fluid nodal pressures are 
used as illustrated below to compute acoustic velocities and intensities. 


Acoustic Intensity 

Intensity is defined as the time averaged product of a pressure with the in-phase component of particle 
velocity. For time-harmonic analysis, where complex numbers are used, this calculation may be visualized as 
taking the dot product of the pressure and velocity phasors. Multiplying one complex number by the in-phase 
part of another complex number is the same operation as multiplying the first number by the complex conjugate 
of the other number and taking the real part of the result: 


where p is pressure, and v is the complex conjugate of velocity. (Often, a factor of 1/2 appears in intensity 
equations. However there is no factor 1/2 in the equations if the assumption is made that pressures and velocities 
are “effective” rms values rather than amplitudes. With this assumption, consistency is maintained, and there is 
no mixing of effective and peak quantities in this formulation.) 


Acoustic Velocities 

The derivation of acoustic velocities for axisymmetric problems is performed for the cosine coefficients 
of the Fourier summation about the axis of rotational symmetry (z), where the r, z, and rotation about 0 DOF 
are active. Only the r and z variations of the scalar pressure field are used to calculate the acoustic velocity vector 
field. The particle velocity in a fluid domain is defined as: 

( 22 ) 


v = 


WQ 


Vp - 


— (*Ef 

o)Q <3r 


*£) 

dz 


where q is the fluid density, w is the circular frequency, iis the square root of -1, and £ are the pressure 


derivatives in the r and z directions, and i and k arc unit vectors in the r and z directions, respectively. A first 
order finite difference approximation of the pressure derivatives at the nodes in an individual 1 RIAAX element 

can be made as shown in Figure 1. The pressure differences between nodes are divided by the distances between 

nodes to approximate the first derivative of pressure in the direction of the two nodes. The approximate pressure 



gradient equations are written for all nodes connected to the node for which the velocity vector is to be calculated. 
In the case above, all elements connected to node 1 must be found, and equations are written for each node 
connected to node 1. No duplicate equations are written when two elements share a common edge. An 
overdetermined system of equations for the pressure derivatives in the r and z directions is the result, with one 
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equation for each node connected to the primary node. The system of equations is 



( 23 ) 


where s* and s* are the r and z components of the unit vector from the primary node to the connecting node 
i, and N is the number of connecting nodes. 

The derivatives are determined approximately using a least squares approach. The particle velocity 
vector is then solved for using Equation 22, and the acoustic intensity vector is given by Equation 21. The 
axisymmetric acoustic intensity field for a complete domain is found by repeating this procedure for each node 
in the domain. In this way, a complete energy flow solution can be derived from nodal pressures and element 
connectivities. 

EXAMPLES 

The axisymmetric model of a submerged, half prolate spheroidal shell is shown in Figure 2. ITie 
structure has a semi axes of 10 m and 5 m. The material is steel with a uniform thickness of 25 mm, a Young s 
modulus (E) of 2.074 Ell N/m 2 , a density (p) of 7860 kg/m 3 , a Poisson’s ratio (v) of 0.3, and a material loss 
factor of 0.0. The frequency range of interest was 100 to 500 Hz. The problem was analyzed for harmonic zero, 
commonly referred to as the “breathing mode” of the domain, implying no variation in the solution field about 
the z axis. 



Figure 2. Axisymmetric Model of a Submerged, Half Prolate-spheroid. 
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As shown in the diagram, the structure is submerged in fluid, which was also modeled and interlaced 
with the structure. The fluid is seawater, with a density (p) ot' 1025 kg/m 3 and a speed of sound (c) of 1500 m/s. 
Two sets of boundary conditions along the left, vertical fluid edge were applied, and arc discussed in more detail 
below. 


A finite element (CONEAX) model of the problem is shown in Figure 3. The structure was modeled 
using axisymmetric conical shell elements. The plate thickness, steel material properties, and the frequency 
range determine an estimated average flexural wavelength [38] of 0.7 m at 500 Hz, which for a mesh requirement 
of about eight elements/wavelength translates to a structural mesh density of 0.087 m/elcment length. 


The outer fluid was modeled using TRI AAX elements. The seawater properties and frequency range 
determine an upper wavelength of 15.0 m and a lower wavelength of 3.0 m. The upper wavelength determines 
the location of the outer boundary of the fluid mesh, shown as A in Figure 2, as one wavelength from the 
structure or 15 0 m. The lower wavelength determines the fluid mesh density, which for a minimum of 8 
elements/wavelength, specifies a fluid mesh size of 0.375 m/element length. The fluid and structural meshes 
are distinct at the fluid-structure boundary, but with coincident nodes. These coincident nodes are coupled by 
area matrices which map fluid pressures to structural forces, as described earlier. Only one fluid DOF was 
assigned to each mesh point. 

The area matrices are input with DMIG cards, which apply area values to the damping, or B2PP matrix. 
A current limitation of NASTRAN is the program performing nonsymmetric system matrix decompositions when 
MPC data is used with DMIG input. Despite the DMIG input being declared symmetric, when the MPC 
equationsare used to obtain the BDD matrix, NASTRAN changes the matrix trailer to nonsymmetric. All matrix 
operations become nonsymmetric as a result, greatly increasing computer time. A sequence ol ALl b 
statements may be used to restore the trailer to its symmetric form. One such sequence (for the 1990 version 
of NASTRAN, Rigid Format 8) is shown below for BDD. 


ALTER 119$ 

DIAGONAL BDD/AVEC/*COLUMN */().$ 
ADD AVEC,/PVEC/(0.0,0.0)$ 

MERGE BDD„„PVEC,/BDDSYM/-1 //6$ 
EQUIV BDDSYM,BDD$ 


Make BDD trailer symmetric 

Vector of ones 

Vector of zeros (P-Vec) 

Dummy merge 

BDD now symmetric 


This alter is inserted before the FRRD module. With the BDD symmetry flag restored, subsequent matrix 
operations will take advantage of the symmetry of the system, reducing the required computer time. his 
NASTRAN bug has been fixed by Gordon Chan of Unisys for the 1992 program release. 

Since the structural elements are about one-quarter the size of the desired fluid elements, some mesh 
transitioning in the fluid meshes was required from the mesh density of the structure (0.087 m/elcment length) 
to the mesh density of the fluid (0.375 m/element length). The fluid clement material properties were assigned 

according to Eq. 19. 


Boundary and Loading Conditions 

Two sets of boundary conditions were applied in the problem. In both cases, a pc impedance, or plane 

wave absorbing boundary condition was applied to the curved outer fluid boundary; and a point forcing function 
was applied at z = 6.0 in the positive r direction along the structure, as shown in Figure 2. The structure was 
constrained in all degrees of freedom at the upper left end. The left vertical fluid boundary was modeled in two 
ways: (1) as a rigid wall, which reflects incident waves, and (2) as a pc impedance absorber. DMIG statements 

inputting -A/pc values for all boundary points simulated the pc plane wave absorber. 


Results 

The problem may be accurately solved for any set of frequencies front 1(X) to 500 I Iz. Detailed results 
are presented here for the 100 Hz case for Cases 1 and 2. The numerical output is examined several ways; (tower 
input and power radiation calculations, structural displacements, acoustic pressure contour plots, and acoustic 
intensity vector plots. This set of output provides a nearly complete solution to the structural-acoustic problem. 
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Figure 3. Finite Element Mesh. 

Since both materials are lossless, all power entering the structure must eventually be radiated into the 
fluid. This means the power input by the forcing function must equal the power exiting the system through the 
plane-wave absorbing outer boundaries. Power input is defined as 

P = F v- ( 24 ) 

1 in 1 in v lih 

where F m is the radial complex input force (in this case unity), and v'„ is the complex conjugate ol the 
corresponding velocity at the force point. The power radiated through the absorbing boundary may be found 
by integrating the calculated acoustic intensity field normal to the boundary: 

P„u, - f V • n dA, ( 25 > 


where C denotes the boundary contour, n is the outward normal vector, and dA is the incremental area. This 
integration may be converted to a summation over all nodes on the boundary, where 

p,, (1 - Y i n, A, (26) 

i - 1 

In this case, n, and A, arc calculated for each node on the boundary, and used with the intensity vector at each 

node to calculate ihe power leaving the system. The summation overall boundary nodes gives the total radiated 
power. 
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Radiated power is commonly represented by a spherical pressure field with a 
yard If the outer boundary is in the far-field, the pressure at one yard is 

' fTTc p~ (27) 

p.yUdB) = 20 log y^r+ m 

where 120 dB is added to calculate the pressure relative to 1 tiPa. 

Thblel lists the power quasi, ties defined above^ and the relate Sefcat 

radiated power from the input power and normatemg toUi S „ as „f a plane-wave 
is probably due to the curved outer boundary. a PF«^ waves is perfectly plane, small reflections at 


Case 

Power Input (W) 

Power Radiated (W) 

%Error 

Pressure@l yard 
(dB re: 1 pPa) 

i 

| 1.499E-8 

1.519E-8 

-1.33 

91.9 

2 

8.988E-8 

8.798E-8 

2.11 

99.5 


Thblel. Power Results for 100 Hz. 

the load. 

Acoustic intensity vector plots are superimposed °"““" a 'Xhe Md^tL^maure 8 “Sure 7 

Figure 5 shows the e "tiie M^or Case l and “ “a Smmon pressure scale is used for all plots, with 

is a plot representing Case 2, and Figure o is a c vector lengths are 

the letters on the contours corresponding to the J^uni eveb (d y e^ ™”“ OT g et ,, e 1/r 

ThP nlntc for Case 1 (Figure 5 and Figure 6) show the highest levels of far-field pressures to be in the 
r and z dhertfonsfuddp^ues between contound and ^urj® ^amt^ng^^e Ibr-fiew^wi^y ^^lORShou^a 

paths in the r and z directions. The ngic 1 boun f d intensity vectors shows an energy path that 

£SbSSS^5=33£E£SS 

^SSSSSSwSsssssfflsssK 

of power though. 
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vectors along the pc boundary now have a z-direcled component, implying power exiting the system through the 
boundary. The overall radiated power is higher for Case 2, though, due to an increase in the radiated energy 

along the midsection of the outer boundary. The near-field intensity plot reveals that the circulating power nows 

between structure and fluid to the left of the load point have disappeared. The energy now flows along the 
fluid-structure boundary and radiates outward at the top of the structure. To the right of the load however, power 
now circulates. The dominant path is from the load point to the right; into and out of the structure; then along 
the structure until r = 0, where some power reenters the structure, and the rest radiates outward along the z axis. 


NASTRAN DIFFICULTIES ENCOUNTERED 

Two important limitations of NASTRAN became apparant during this study: the BDD damping matrix's 
being specified as unsymmetric regardless of the symmetric nature of fluid structure interaction and absorbing 
boundary data input by DMIG cards; and the formulation of stiffness coefficients for CONEAX elements. The 
BDD matrix trailer maybe restored to symmetric using the ALTER statements outlined in the Example section. 
The difficulties with the CONEAX formulation are not easily fixed however. 

Stiffnesses for CONEAX elements are computed analytically by NASTRAN, and involve the inverse 
of Ar/Al .where Ar is the difference in radii and A1 is the total distance between the grids defining an element. 
For perfectly cylindrical shell elements with no variation in radii (Ar/Al = 0), a different formulation is used 
to avoid a floating point error caused by a division by zero. However, no provision is made for small relative 
variation in radii ( Ar/Al ~ 0), and for a small range of elements the analytical computation is corrupted when 
computer precision limits are reached. Sometimes the error is so drastic that negative values are obtained for 
self term (diagonal) stiffnesses. The negative stiffnesses are reported to the user when NASTRAN checks the 
system matrices for singularities. However, sometimes the error may be drastic in the positive sense, i.e., 
stiffnesses orders of magnitude too large. No error would be reported to the user, and the final solution would 
be incorrect. 


CONEAX stiffness errors were encountered for the example described here at the upper left end of 
the structure, where the elements become nearly cylindrical. In this case, the stiffnesses ol several oi the near 
cylindrical elements were output and analyzed for accuracy. The two end elements were found to have large 
errors in stiffness. To solve the problem, the radii of the element grid points were set equal, and the end of the 
structure was approximated as purely cylindrical. 

A possible programming solution to the sensitivity of CONEAX stiffnesses to small relative differences 
in radii is to approximate nearly cylindrical regions as cylindrical. For example, if for a given element Ar/Al 

is below some specified tolerance e, the second grid radius is set equal to the first grid radius. The resulting model 
would be a stepwise approximation of the nearly cylindrical region. The chief problem is how to choose e. Studies 
would have to be performed on ranges of nearly cylindrical elements using different levels of computer precision 
to determine the accuracy limits on the analytical stiffness computation method. 


CONCLUSIONS 

The combination of structural displacement plots, pressure contours, and acoustic intensity vector fields 
all serve to reveal the complete state of a structural-acoustic problem. However, one component of the response 
is missing: the energy flow within the structure. The circulating energy along the structural-acoustic boundary 
indicated by the intensity plots show power flowing through the structure. A formulation similar to that for 
acoustic intensity can be performed for the structure; however more than one wave type must be considered. 
For the axisymmetric shells of revolution (CONEAX) used here, for example, both flexural (composed of both 
shear and moment waves) and longitudinal waveforms may transport energy through structures. Methods have 
been developed for general three-dimensional structural models of beams (BAR) and plates (QUAD4) [39], but 
have not yet been extended to axisymmetric problems. This additional analysis tool will help improve 
considerably the understanding of structural-acoustic/ frequency response problems. 
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Undisplaced Shape 

Displaced Shape (rigid Boundary Condition) 
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Figure 4. Structural Displacements for both Boundary Condition Cases (phase = 292.5°) 




l igme 5. Acoustic Intensities 
I Hies Denote Constant Pressure 


and Pressures for Rigid Vertical Boundary; 
Contours; Vectors Denote Acoustic Intensities. 
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Fieure 6. Close-up of Acoustic Intensities and Pressures for Rigid Vertical Boundary; 
Lines Denote Constant Pressure Contours; Vectors Denote Acoustic Intensities. 
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Figure 7. Acoustic Intensities and Pressures for Absorbing Vertical Boundary, 
Lines Denote Constant Pressure Contours; Vectors Denote Acoustic Intensities 
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